What is LiDAR?

LiDAR (Light detection and ranging) is a data acquisition technology that uses light emitted from a laser pulse (up to 400000 pulses of light per second) to measure the travel time of the pulse from the source to the target and back (returns). Using this approach the instrument collects three dimensional data in large volumes, high density and at unprecedented accuracy (Figure 1). The resolution of the data collected depends on the pulse frequency, that is, resolution increases with increasing pulse frequency.

What types of variables can be extracted from LiDAR data?

LiDAR data is provided in LAS format and classified into ground and non-ground (vegetation) returns. This type of classification is usually performed by the data provider.

  1. Canopy density map -. Forest canopy density is the ratio of the points classified as vegetation to the point classified as ground. This variable, together with canopy height, are commonly used for biomass estimation and the assessment of vegetation coverage.

  2. Digital Terrain Model (DTM) -. The DTM represents the bare ground, it can be calculated from the average value of points that have been classified as ground in each cell.

  3. Canopy height metrics -. The forest canopy height is calculated by substracting the DTM from the first return. Other canopy height metrics such as maximum, minimum, mean, standard deviation, variance and percentiles can also be calculated.

    LiDAR data processing

Setting the R environment

The package rgrass7 provides the interface between GRASS GIS 7 and the R software that allows the use of all the GRASS GIS commands from the R command line. To load the rgrass7 package into your R session, it is necessary to first install the package using the command install.package(“rgrass7”) and then:

library(rgrass7)
## Loading required package: XML
## GRASS GIS interface loaded with GRASS version: GRASS 7.0.4svn (2016)
## and location: DIARS

It is recommended to first create a folder where the datasets (LAS files, download a sample dataset from HERE) are located and the outputs and results will be saved. To set the working directory in R to the folder created:

setwd("/home/garzonc/Desktop/DIARS/")

Now you need to create a directory that will contain the GRASS GIS database for the working session:

# dir.create('grassdata')

To initiate GRASS into your R session you need to set the gisBase argument to GRASS binaries and libraries directory path, this should be something like C:\Program Files\GRASS GIS 7.00beta4 if your are using Windows (check HERE for other operating systems). This will create a new GRASS location called Sylt_island inside the GRASS GISDBASE.

myGRASS <- "/usr/local/grass-7.0.4svn"
myGISDbase <- "/home/garzonc/Documents/grassData"
myLocation <- "DIARS"
myMapset <- "PERMANENT"

initGRASS(myGRASS, home = tempdir(), gisDbase = myGISDbase, 
    location = myLocation, mapset = myMapset, 
    override = TRUE)
## gisdbase    /home/garzonc/Documents/grassData 
## location    DIARS 
## mapset      PERMANENT 
## rows        1440 
## columns     1440 
## north       6075220 
## south       6074500 
## west        453999.9 
## east        454719.9 
## nsres       0.5 
## ewres       0.5 
## projection  +proj=utm +no_defs +zone=32
## +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000
## +to_meter=1

And now set the projection system, for this example we need WGS84 UTM zone 32 (epgs: 32632). You can search for the epgs code of your projection system here:

execGRASS("g.proj", flags = c("c"), parameters = list(epsg = 32632))
## Projection information updated

To set the geographic region in GRASS including the extent of the bounding box (i.e. the region where all the calculations will be performed) and the resolution:

execGRASS("g.region", parameters = list(n = "6075220.40", 
    s = "6074500.40", e = "454719.89", w = "453999.89", 
    res = "0.5"))

To make sure that the GRASS location is correct you can check the location metadata:

gmeta()
## gisdbase    /home/garzonc/Documents/grassData 
## location    DIARS 
## mapset      PERMANENT 
## rows        1440 
## columns     1440 
## north       6075220 
## south       6074500 
## west        453999.9 
## east        454719.9 
## nsres       0.5 
## ewres       0.5 
## projection  +proj=utm +no_defs +zone=32
## +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000
## +to_meter=1

Importing the LiDAR data into R via GRASS GIS

Once the GRASS in R environment is set you can import the LiDAR data. The following for loop will create 4 vector point files from the binary LAS files using the v.in.lidar GRASS module. Note that this will take a couple of hours to run depending on your computer configuration because this files contain altogether 16568266 points.

# for (i in
# c('GLOBAL-H_454000_6074500.las','GLOBAL-H_454000_6075000.las','GLOBAL-H_454500_6074500.las','GLOBAL-H_454500_6075000.las'))
# {execGRASS('v.in.lidar',
# flags=c('r','o','overwrite','quiet'),
# parameters=list(input=i,output=substr(i,8,23)))}

Now combine the 4 imported vector map layers into one single vector map layer (this step is also a time consuming task):

# execGRASS('v.patch', flags=c('e',
# 'overwrite','quiet'),
# parameters=list(input=c('H_454000_6074500','H_454000_6075000','H_454500_6074500','H_454500_6075000'),output='las_file'))

Select the points classified as “vegetation” and create a new vector map:

# execGRASS('v.extract',flags=c('overwrite','quiet'),
# parameters=list(input='las_file',where='class
# IN ('Low Vegetation','High
# Vegetation','Medium
# Vegetation')',output='las_file_vegetation'))

Select the points classified as “ground” and create a new vector map:

# execGRASS('v.extract',flags=c('overwrite','quiet'),
# parameters=list(input='las_file',where='class
# IN
# ('Ground')',output='las_file_ground'))

Select the “non-vegetation” points and create a new vector map:

# execGRASS('v.extract',flags=c('overwrite','quiet'),
# parameters=list(input='las_file',where='class
# IN ('Unclassified','Ground','Building'
# ,'Water')',output='las_file_non_vegetation'))

Convert each of the 3 GRASS binary vector map containing the vegetation, ground and non-vegetation points into a GRASS ASCII vector map:

# execGRASS('v.out.ascii',flags='overwrite',parameters=list(input='las_file_vegetation',output='las_file_vegetation',format='point'))
# execGRASS('v.out.ascii',flags='overwrite',parameters=list(input='las_file_non_vegetation',output='las_file_non_vegetation',format='point'))
# execGRASS('v.out.ascii',flags='overwrite',parameters=list(input='las_file_ground',output='las_file_ground',format='point'))

To visualise a subset of the point cloud using the RGL package you need to first install the rgl package using the install.package("rgl") command:

library(rgl)
ground <- read.table("las_file_ground", sep = "|", 
    dec = ".")
ground <- ground[which(ground[, 1] < 454100 & 
    ground[, 2] > 6074800 & ground[, 2] < 
    6074900), ]
plot3d(ground[, 1], ground[, 2], ground[, 
    3], col = "coral", aspect = c(1, 1, 0.1), 
    xlab = "", ylab = "", zlab = "", box = F, 
    axes = F)
vegetation <- read.table("las_file_vegetation", 
    sep = "|", dec = ".")
vegetation <- vegetation[which(vegetation[, 
    1] < 454100 & vegetation[, 2] > 6074800 & 
    vegetation[, 2] < 6074900), ]
plot3d(vegetation[, 1], vegetation[, 2], 
    vegetation[, 3], col = "green", aspect = c(1, 
        1, 0.1), add = T)

You must enable Javascript to view this page properly.

Creating a canopy density map

Create a raster containing the number of points classified as vegetation in each cell:

execGRASS("r.in.xyz", flags = "overwrite", 
    parameters = list(input = "las_file_vegetation", 
        output = "raster_vegetation", method = "n", 
        type = "CELL"))
## Reading input data...
##   18%  37%  55%  74% 100%
## Writing to output raster map...
##    0%   6%  12%  18%  24%  30%  36%  42%  48%  54%  60%  66%  72%  78%  84%  90%  96% 100%
## r.in.xyz complete. 497051 points found in region.

Create a raster containing the number of non-vegetation points in each cell:

execGRASS("r.in.xyz", flags = "overwrite", 
    parameters = list(input = "las_file_non_vegetation", 
        output = "raster_non_vegetation", 
        method = "n", type = "CELL"))
## Reading input data...
##    3%   7%  11%  15%  19%  23%  27%  31%  35%  39%  43%  47%  51%  55%  59%  63%  67%  71%  75%  79%  83% 100%
## Writing to output raster map...
##    0%   6%  12%  18%  24%  30%  36%  42%  48%  54%  60%  66%  72%  78%  84%  90%  96% 100%
## r.in.xyz complete. 15511478 points found in region.

Convert any resulting NoData cells to 0 so that the subsequent operations treat NoData cells as 0:

execGRASS("r.null", parameters = list(map = "raster_vegetation", 
    null = 0))
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%
execGRASS("r.null", parameters = list(map = "raster_non_vegetation", 
    null = 0))
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%

Canopy density can be calculated as the percentage of total returns that were vegetation returns. This gives us the ratio from 0 to 1 where 0 represents no vegetation and 1 very dense vegetation cover.

execGRASS("r.mapcalc", flags = "overwrite", 
    parameters = list(expression = "Canopy_density=double(raster_vegetation)/(double(raster_vegetation)+double(raster_non_vegetation))"))

To import the density canopy map from GRASS to R and plot the result:

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## Creating BIL support files...
## Exporting raster as floating values (bytes=8)
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
spplot(density.map)

Creating a Digital Terrain Model (DTM)

A DTM can be created by calulating the average value of points that represent the bare earth surface in each cell (cf. classified as “ground”).

execGRASS("r.in.xyz", flags = "overwrite", 
    parameters = list(input = "las_file_ground", 
        output = "DTM", method = "mean", 
        type = "CELL"))
## Reading input data...
##    3%   7%  11%  15%  19%  23%  27%  31%  35%  39%  43%  47%  51%  55%  59%  63%  67%  71%  75%  79%  83% 100%
## Writing to output raster map...
##    0%   6%  12%  18%  24%  30%  36%  42%  48%  54%  60%  66%  72%  78%  84%  90%  96% 100%
## r.in.xyz complete. 15511478 points found in region.

To import the DTM from GRASS to R and plot the result:

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## Creating BIL support files...
## Exporting raster as integer values (bytes=4)
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
spplot(DTM)

Note that some cells have NA values, in this case it will be better to use the v.surf.idw GRASS module to interpolate elevation in these pixels using the inverse distance squared weighting (IDW).

execGRASS("v.surf.idw", flags = "overwrite", 
    parameters = list(input = "las_file_ground", 
        output = "DTM", column = "z_coord"))
## 15513638 points loaded
## Interpolating raster map <DTM> (1440 rows, 1440 columns)...
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%
## v.surf.idw complete.
DTM <- readRAST("DTM")
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## Creating BIL support files...
## Exporting raster as floating values (bytes=8)
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
spplot(DTM)

Creating rasters with canopy height metrics

To calculate canopy height metrics raster files (maximum, minimum, mean, standard deviation,variance, 5th and 95th percentile):

execGRASS("r.in.xyz", flags = "overwrite", 
    parameters = list(input = "las_file_vegetation", 
        output = "Max_height", method = "max"))
## Reading input data...
##   18%  37%  55%  74% 100%
## Writing to output raster map...
##    0%   6%  12%  18%  24%  30%  36%  42%  48%  54%  60%  66%  72%  78%  84%  90%  96% 100%
## r.in.xyz complete. 497051 points found in region.
execGRASS("r.in.xyz", flags = "overwrite", 
    parameters = list(input = "las_file_vegetation", 
        output = "Min_height", method = "min"))
## Reading input data...
##   18%  37%  55%  74% 100%
## Writing to output raster map...
##    0%   6%  12%  18%  24%  30%  36%  42%  48%  54%  60%  66%  72%  78%  84%  90%  96% 100%
## r.in.xyz complete. 497051 points found in region.
execGRASS("r.in.xyz", flags = "overwrite", 
    parameters = list(input = "las_file_vegetation", 
        output = "Mean_height", method = "mean"))
## Reading input data...
##   18%  37%  55%  74% 100%
## Writing to output raster map...
##    0%   6%  12%  18%  24%  30%  36%  42%  48%  54%  60%  66%  72%  78%  84%  90%  96% 100%
## r.in.xyz complete. 497051 points found in region.
execGRASS("r.in.xyz", flags = "overwrite", 
    parameters = list(input = "las_file_vegetation", 
        output = "Stddev_height", method = "stddev"))
## Reading input data...
##   18%  37%  55%  74% 100%
## Writing to output raster map...
##    0%   6%  12%  18%  24%  30%  36%  42%  48%  54%  60%  66%  72%  78%  84%  90%  96% 100%
## r.in.xyz complete. 497051 points found in region.
execGRASS("r.in.xyz", flags = "overwrite", 
    parameters = list(input = "las_file_vegetation", 
        output = "Variance_height", method = "variance"))
## Reading input data...
##   18%  37%  55%  74% 100%
## Writing to output raster map...
##    0%   6%  12%  18%  24%  30%  36%  42%  48%  54%  60%  66%  72%  78%  84%  90%  96% 100%
## r.in.xyz complete. 497051 points found in region.
execGRASS("r.in.xyz", flags = "overwrite", 
    parameters = list(input = "las_file_vegetation", 
        output = "Percentile95_height", method = "percentile", 
        pth = 95))
## Reading input data...
##   18%  37%  55%  74% 100%
## Writing to output raster map...
##    0%   6%  12%  18%  24%  30%  36%  42%  48%  54%  60%  66%  72%  78%  84%  90%  96% 100%
## r.in.xyz complete. 497051 points found in region.
execGRASS("r.in.xyz", flags = "overwrite", 
    parameters = list(input = "las_file_vegetation", 
        output = "Percentile5_height", method = "percentile", 
        pth = 5))
## Reading input data...
##   18%  37%  55%  74% 100%
## Writing to output raster map...
##    0%   6%  12%  18%  24%  30%  36%  42%  48%  54%  60%  66%  72%  78%  84%  90%  96% 100%
## r.in.xyz complete. 497051 points found in region.

To substract the DEM from each of the canopy height rasters to produce a canopy height layers normalized for topography:

exp1 <- "Max_height=if(double(Max_height)<double(DTM),0,double(Max_height)-double(DTM))"
execGRASS("r.mapcalc", flags = "overwrite", 
    parameters = list(expression = exp1))
exp2 <- "Min_height=if(double(Min_height)<double(DTM),0,double(Min_height)-double(DTM))"
execGRASS("r.mapcalc", flags = "overwrite", 
    parameters = list(expression = exp2))
exp3 <- "Mean_height=if(double(Mean_height)<double(DTM),0,double(Mean_height)-double(DTM))"
execGRASS("r.mapcalc", flags = "overwrite", 
    parameters = list(expression = exp3))
exp4 <- "Stddev_height=if(double(Stddev_height)<double(DTM),0,double(Stddev_height)-double(DTM))"
execGRASS("r.mapcalc", flags = "overwrite", 
    parameters = list(expression = exp4))
exp5 <- "Variance_height=if(double(Variance_height)<double(DTM),0,double(Variance_height)-double(DTM))"
execGRASS("r.mapcalc", flags = "overwrite", 
    parameters = list(expression = exp5))
exp6 <- "Percentile95_height=if(double(Percentile95_height)<double(DTM),0,double(Percentile95_height)-double(DTM))"
execGRASS("r.mapcalc", flags = "overwrite", 
    parameters = list(expression = exp6))
exp7 <- "Percentile5_height=if(double(Percentile5_height)<double(DTM),0,double(Percentile5_height)-double(DTM))"
execGRASS("r.mapcalc", flags = "overwrite", 
    parameters = list(expression = exp7))

To convert any resulting NoData cells to 0:

execGRASS("r.null", parameters = list(map = "Max_height", 
    null = 0))
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%
execGRASS("r.null", parameters = list(map = "Min_height", 
    null = 0))
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%
execGRASS("r.null", parameters = list(map = "Mean_height", 
    null = 0))
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%
execGRASS("r.null", parameters = list(map = "Stddev_height", 
    null = 0))
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%
execGRASS("r.null", parameters = list(map = "Variance_height", 
    null = 0))
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%
execGRASS("r.null", parameters = list(map = "Percentile95_height", 
    null = 0))
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%
execGRASS("r.null", parameters = list(map = "Percentile5_height", 
    null = 0))
##    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%

To import the maximum canopy height map from GRASS to R and plot the result:

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## Creating BIL support files...
## Exporting raster as floating values (bytes=8)
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
spplot(Max_height)

Back to the top

Authors: Tarek Hattab, Carol Garzon-Lopez, Duccio Rocchini & Jonathan Lenoir

Date: 2016-04-15